1.Import packages

If there are errors take place, you can run install.packages({missing package name}) to install packages.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✔ ggplot2 3.2.1     ✔ purrr   0.3.3
## ✔ tibble  3.0.3     ✔ dplyr   1.0.2
## ✔ tidyr   1.1.2     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following object is masked from 'package:purrr':
## 
##     transpose
library(sp)
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
library(tm)
## Loading required package: NLP
## 
## Attaching package: 'NLP'
## The following object is masked from 'package:ggplot2':
## 
##     annotate
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Loading required package: rpart
## Registered S3 method overwritten by 'spatstat':
##   method     from
##   print.boxx cli
## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: spatstat version 1.64-1 is out of date by more than 8 months; we recommend upgrading to the latest version.
## 
## Attaching package: 'spatstat'
## The following object is masked from 'package:data.table':
## 
##     shift
library(here)
## here() starts at /Users/fangzeqiang/Desktop/SDSV/Final/CASA0005_Final_Assessment
library(sp)
library(rgeos)
## rgeos version: 0.5-5, (SVN revision 640)
##  GEOS runtime version: 3.7.2-CAPI-1.11.2 
##  Linking to sp version: 1.4-2 
##  Polygon checking: TRUE
library(maptools)
## Checking rgeos availability: TRUE
## 
## Attaching package: 'maptools'
## The following object is masked from 'package:table1':
## 
##     label
library(tmap)
library(sf)
library(geojson)
## 
## Attaching package: 'geojson'
## The following object is masked from 'package:graphics':
## 
##     polygon
library(geojsonio)
## 
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
## 
##     pretty
library(tmaptools)
library(RColorBrewer)
library(spdep) 
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

2. The basic settings

2.1 Set the path of your project.

The path below is mine, you should set your own work path

setwd("/Users/fangzeqiang/Desktop/SDSV/Final/CASA0005_Final_Assessment/")

2.2 Import the shape file

What you should keep in mind is that this shape file should be run in the complete ESRI dir

# London_Borough=st_read("F:/spat/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp")
London_Borough = st_read("dataset/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp")
## Reading layer `London_Borough_Excluding_MHW' from data source `/Users/fangzeqiang/Desktop/SDSV/Final/CASA0005_Final_Assessment/Dataset/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp' using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## CRS:            27700
# plot the map
plot(st_geometry(London_Borough))

2.3 import the London national chargepoint register dataset

df = fread("Dataset/London_NCR_GSS_Added.csv")

3. The distribution of samples

#3.1 select the years
df$year = year(df$dateCreated)
table(df$year)
## 
## 2017 2018 2019 2020 2021 
##    1    4  696  362    3
#select 2019 & 2020 
df = df[year==2019|year==2020,]
table(df$year)
## 
## 2019 2020 
##  696  362
#3.2 show the table of the distribution of two years samples
table1(~county|factor(year),data=df)
2019
(N=696)		 2020
(N=362)		 Overall
(N=1058)
county
London 150 (21.6%) 150 (41.4%) 300 (28.4%)
London Borough of Camden 87 (12.5%) 179 (49.4%) 266 (25.1%)
London Borough of Ealing 65 (9.3%) 1 (0.3%) 66 (6.2%)
London Borough of Greenwich 2 (0.3%) 0 (0%) 2 (0.2%)
London Borough of Hackney 62 (8.9%) 0 (0%) 62 (5.9%)
London Borough of Hammersmith and Fulham 2 (0.3%) 3 (0.8%) 5 (0.5%)
London Borough of Hounslow 78 (11.2%) 0 (0%) 78 (7.4%)
London Borough of Islington 3 (0.4%) 15 (4.1%) 18 (1.7%)
London Borough of Lambeth 118 (17.0%) 2 (0.6%) 120 (11.3%)
London Borough of Richmond upon Thames 8 (1.1%) 9 (2.5%) 17 (1.6%)
London Borough of Southwark 1 (0.1%) 0 (0%) 1 (0.1%)
London Borough Of Southwark 57 (8.2%) 0 (0%) 57 (5.4%)
London Borough of Waltham Forest 60 (8.6%) 2 (0.6%) 62 (5.9%)
London Borough of Wandsworth 3 (0.4%) 1 (0.3%) 4 (0.4%) 
head(df)
##                      chargeDeviceID latitude longitude      town county
## 1: b353888fbc3ee702b41d30669a23e12d 51.54491  -0.00876 Stratford London
## 2: f53d6d41d4ea910e11d4ea914d58b803 51.54491  -0.00876 Stratford London
## 3: 61ef7c5a4cdc44aa15d46293f1f185b8 51.54491  -0.00876 Stratford London
## 4: 3ae50acb1eb21a62b4f1d6620f155c10 51.54491  -0.00876 Stratford London
## 5: 3c348a8f6ea52056921b931828359346 51.54491  -0.00876 Stratford London
## 6: 6e785fdfb123ff8b192f26cf595f74b1 51.54491  -0.00876 Stratford London
##    postcode chargeDeviceStatus         dateCreated         dateUpdated
## 1:  E20 1YY         In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 2:  E20 1YY         In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 3:  E20 1YY         In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 4:  E20 1YY         In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 5:  E20 1YY         In service 2020-01-06 15:36:08 2020-01-06 15:36:08
## 6:  E20 1YY         In service 2020-01-06 15:36:08 2020-01-06 15:36:08
##    lastUpdated    locationType  GSS_CODE year
## 1:             Public car park E09000025 2020
## 2:             Public car park E09000025 2020
## 3:             Public car park E09000025 2020
## 4:             Public car park E09000025 2020
## 5:             Public car park E09000025 2020
## 6:             Public car park E09000025 2020

4. Visualisation and comparison of the density of London EV charge point between two years

4.1 Data cleaning for mapping

  • I divided the dataset into ones in two years, and merged these two processed dataset with geometric data.
  • Then we calculate the density for two years.
  • Finally we can visualise and map the density
#5.1 process the data to divide them by years(2019 & 2020)
df1<-subset(df,year==2019) #2019 year
df2<-subset(df,year==2020) #2020 year

# EV charge points created in 2019
# sdf1<-merge(London_Borough,df1,by="GSS_CODE")
sdf1<-merge(London_Borough,df1,by="GSS_CODE",all = TRUE)
sdf1<-sdf1[,c("GSS_CODE","geometry","longitude","latitude")]

# EV charge points created in 2020
# sdf2<-merge(London_Borough,df2,by="GSS_CODE")
sdf2<-merge(London_Borough,df2,by="GSS_CODE",all = TRUE)
sdf2<-sdf2[,c("GSS_CODE","geometry","longitude","latitude")]

### 4.2 The density of data in 2019

# Data preparation
nsdf1 = sdf1%>%
  add_count(GSS_CODE)%>%
  mutate(area=st_area(.))%>%
  # Use dplyr::mutate to calculate the density of the charge point for each borough
  mutate(density=n*1000*1000/area)
  # because the st_area default unit is square metre

# select the following variables---"density","GSS_CODE","n"(the count of GSS_CODE)
nsdf1 = dplyr::select(nsdf1,density,GSS_CODE, n)

nsdf1 = nsdf1%>%                    
  group_by(GSS_CODE)%>%         
  summarise(density =first(density),GSS_CODE=first(GSS_CODE))
## `summarise()` ungrouping output (override with `.groups` argument)
tmap_mode("plot")
## tmap mode set to plotting
# plot the figure2: The distribution of the density of the London charge points in 2019
tm_shape( nsdf1) +
  tm_compass( north = 0,
              type = "4star",
              text.size = 0.8,
              size = 2.5,
              show.labels = 1,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              position = c("left","top"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA,
              fontsize = 1.5) +
  tm_polygons("density",
              style="jenks",
              palette="RdPu",
              midpoint=NA,
              popup.vars=c("GSS_CODE", "density"),
              title="Density per square kilometre (2019)"
              )
## Warning: The argument fontsize of tm_compass is deprecated. It has been renamed
## to text.size

Figure the London boroughs map with their names

### 4.3 The density of data in 2020

nsdf2 = sdf2%>%
  add_count(GSS_CODE)%>%
  mutate(area=st_area(.))%>%
  mutate(units::set_units(area,km^2))%>%
  #then density of the points per ward
  mutate(density=n*1000*1000/area)

nsdf2 = dplyr::select(nsdf2,density,GSS_CODE, n)

nsdf2 = nsdf2%>%                    
  group_by(GSS_CODE)%>%         
  summarise(density =first(density),GSS_CODE=first(GSS_CODE))
## `summarise()` ungrouping output (override with `.groups` argument)
tmap_mode("plot")
## tmap mode set to plotting
tm_shape( nsdf2) +
  tm_compass( north = 0,
              type = "4star",
              text.size = 0.8,
              size = 2.5,
              show.labels = 1,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              position = c("left","top"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA,
              fontsize = 1.5) +
  tm_polygons("density",
              style="jenks",
              palette="PuOr",
              midpoint=NA,
              popup.vars=c("GSS_CODE", "density"),
              title="Density per square kilometre (2020)")
## Warning: The argument fontsize of tm_compass is deprecated. It has been renamed
## to text.size

5. Analysing Spatial Autocorrelation with Moran’s I

Since the sample in 2019 & 2020 is too small to run a good result, I analysis the autocorrelation based on the samples which contain these two yearsl

###5.1 generate the data for analysis

sdf = merge(London_Borough,df,by="GSS_CODE",all = TRUE)
sdf = sdf[,c("GSS_CODE","geometry","longitude","latitude")]
nsdf = sdf%>%
  add_count(GSS_CODE)%>%
  mutate(area=st_area(.))%>%
  #then density of the points per ward
  mutate(density=n*1000*1000/area)

nsdf = dplyr::select(nsdf,density,GSS_CODE, n)

nsdf = nsdf%>%                    
  group_by(GSS_CODE)%>%         
  summarise(density = first(density), GSS_CODE = first(GSS_CODE))
## `summarise()` ungrouping output (override with `.groups` argument)

###5.2 First plot the centroids of all boroughs in London

coordsW = nsdf%>%
   st_centroid()%>%
   st_geometry()
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x
 plot(coordsW,axes=TRUE)

5.3 create a neighbours list

 LWard_nb <- nsdf%>%poly2nb(.,queen=T)

plot the neighbours list we creat

plot(nsdf$geometry)
plot(LWard_nb, st_geometry(coordsW), col="blue", add = T)

# plot(LWard_nb, st_geometry(coordsW), col="red")
# plot(nsdf$geometry, add=T)
 # plot(LWard_nb, st_geometry(coordsW), col="red")
 # jpeg(file="saving_plot1.jpeg")
 # dev.off()

### 5.3 create a spatial weights object from these weights

Lward.lw = nb2listw(LWard_nb, style="C")
head(Lward.lw$neighbours)
## [[1]]
## [1]  7 12 19 30 33
## 
## [[2]]
## [1] 16 25 26
## 
## [[3]]
## [1]  5  7 10 14 15
## 
## [[4]]
## [1]  6 11
## 
## [[5]]
## [1]  3  7  9 13 15 20 33
## 
## [[6]]
## [1]  4  8 11 22 23 28

### 5.4 Calculate the Global Moran’I Index

  • Conduct the global Moran index
 I_LWard_Global_Density = nsdf %>%
   pull(density) %>%
   as.vector()%>%
   moran.test(.,Lward.lw)
names(I_LWard_Global_Density)
## [1] "statistic"   "p.value"     "estimate"    "alternative" "method"     
## [6] "data.name"
head(I_LWard_Global_Density)
## $statistic
## Moran I statistic standard deviate 
##                         0.03513093 
## 
## $p.value
## [1] 0.4859877
## 
## $estimate
## Moran I statistic       Expectation          Variance 
##      -0.028283348      -0.031250000       0.007131055 
## 
## $alternative
## [1] "greater"
## 
## $method
## [1] "Moran I test under randomisation"
## 
## $data.name
## [1] ".  \nweights: Lward.lw    \n"
  • Conduct the Local Moran test
 I_LWard_Local_Density = nsdf %>%
   pull(density) %>%
   as.vector()%>%
   localmoran(., Lward.lw)%>%
   as_tibble()
  • Merge the moran test result with the geometric dataset
nsdf<-nsdf%>%
   mutate(density_I = as.numeric(I_LWard_Local_Density$Ii))%>%
   mutate(density_Iz =as.numeric(I_LWard_Local_Density$Z.Ii))

5.5 Interactive visulisation of the distribution of the local Moran result

tmap_mode("view")
## tmap mode set to interactive viewing
#set the group and colour
summary(nsdf$density_Iz)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.813183 -0.007425  0.133178  0.046786  0.449411  0.975067
breaks1 = seq(-3,1,0.5) 
# breaks2 = c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)
# Depends on the max and min value in Moran's I
MoranColours = rev(brewer.pal(8, "RdGy"))

# Plot the map
tm_shape(nsdf) +
  tm_polygons("density_Iz",
              style="fixed",
              breaks=breaks1,
              palette=MoranColours,
              midpoint=NA,
              title="Local Moran's I,EV charge points in London")

### 5.6 GI score

Gi_LWard_Local_Density = nsdf %>%
  pull(density) %>%
  as.vector()%>%
  localG(., Lward.lw)

head(Gi_LWard_Local_Density)
## [1]  2.3573180 -0.7667792  1.8494548 -0.6694542  1.1946019  0.6283400
Gi_nsdf = nsdf %>%
  mutate(density_G = as.numeric(Gi_LWard_Local_Density))
GIColours = rev(brewer.pal(8, "RdBu"))

#now plot on an interactive map
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(Gi_nsdf) +
    tm_polygons("density_G",
        style="fixed",
        breaks=breaks1,
        palette=GIColours,
        midpoint=NA,
        title="Gi*, EV charge points in London")
## Warning: Values have found that are higher than the highest break

The END.

Thanks for watching!